      SUBROUTINE RDMDL(ISET,JSET,IPMLS)
C
C PURPOSE -   TO READ AN EARTH MODEL.
C
C             THIS PROGRAM READS AN EARTH MODEL, DIVIDES ONE STEP INTO
C     NDIV STEPS UP TO DMAX DEPTH, AND CONVERT VPI AND VSI TO LAMBDA
C     AND MU.
C
C MODEL -     AN EARTH MODEL IS SPECIFIED BY SIX PHYSICAL PARAMETERS,
C     RHOI, VPI, VSI, XII, PHII, AND ETAI AT SUCCESSIVE DEPTHS.
C     WHEN THE MODEL IS ISOTROPIC (ISO=0) XII, PHII, AND ETAI SHOULD NOT
C     BE GIVEN.
C             DEPTH INFORMATION IS GIVEN BY THE ARRAY HI.
C     THE FIRST HI CONTAINS THE DEPTH (USUALLY ZERO) THAT CORRESPONDS TO
C     THE PHYSICAL PARAMETERS GIVEN BY THE FIRST CARD.
C     I-TH ELEMENT OF HI,E.G. HI(I), SPECIFIES THE DISTANCE BETWEEN
C     (I-1)-TH AND I-TH CARD.     IF(HI(I).EQ.0) TWO ORDINATES ARE
C     INTERPRETED AS UPPER(I-1) AND LOWER(I) SIDES OF A DISCONTINUITY.
C
      REAL*8  RAD,TMASS,GS,RHO0,GA
      REAL*8  D,H,RHO,ELAMB,EMU,GR,XI,PHI,ETA
      COMMON/MODEL /NLAY,ISO,LAME,NDIV,DMAX,ILAY,
     1              RAD,TMASS,GS,RHO0,GA,
     2              HI(500),RHOI(500),VPI(500),VSI(500),NAME(20),
     3              XII(500),PHII(500),ETAI(500),
     4              D(1000),H(1000),RHO(1000),ELAMB(1000),EMU(1000),
     5              GR(1000),XI(1000),PHI(1000),ETA(1000)
      REAL*8  DD,HH
C*    DIMENSION  FMT(20)
      CHARACTER*80 FMT
C
C INPUT -     ILAY  =NO. OF ORDINATES(DEPTHS) AT WHICH PARAMETERS ARE
C             GIVEN
C             RAD   =RADIUS OF THE EARTH IN KM
C             TMASS =TOTAL MASS OF THE EARTH IN C.G.S.
C             GS    =SURFACE VALUE OF THE GRAVITY IN C.G.S.
C             RHO0  =DENSITY IN THE INNERMOST HOMOGENENEOUS CORE
C             ISO   =0 WHEN THE MODEL IS ISOTROPIC
C             LAME  =0 WHEN VPI AND VSI ARE VELOCITIES.
C             IF(LAME.NE.0)  THEN VPI AND VSI ARE INTERPRETED AS LAME'S
C             CONSTANTS.
C             HI    =DISTANCES BETWEEN TWO ORDINATES IN KM
C             RHOI  =DENSITIES IN C.G.S.
C             VPI   =COMPRESSIONAL WAVE VELOCITIES IN KM/SEC OR
C                    LAMBDA                        IN (10**10) C.G.S.
C             VSI   =SHEAR WAVE VELOCITIES         IN KM/SEC  OR
C                    MU                            IN (10**10) C.G.S.
C             XII   =(C11-C12)/2C44
C             PHII  =C33/C11
C             ETAI  =C13/(C11-2C44)
C             JSET  =DATA SET REF. NO. FOR OUTPUT
C             ISET  =DATA SET REF. NO. FOR INPUT
C
C OUTPUT -    NLAY  =NO. OF ORDINATES IN DIVIDED MODEL
C             D     =DEPTHS IN KM
C             H     =DISTANCES IN KM     H(1)=D(1),  H(I)=D(I)-D(I-1)
C             RHO   =DENSITIES IN C.G.S.
C             ELAMB =LAMBDA IN (10**10) C.G.S.
C             EMU   =MU     IN (10**10) C.G.S.
C             GR    =GRAVITIES IN (10**5) C.G.S.
C             XI    =XI
C             PHI   =PHI
C             ETA   =ETA
C     IF IPMLS .EQ. 0 PRINT MODEL
C SUBROUTINE REQUIRED - DIVIDE, GRAVTY
C
   
      IF(ISET-7)  100,100,500
C
C CARD INPUT
C
  100 CONTINUE
      READ(5,110)  NAME
     
  110 FORMAT(20A4)
      READ(5,120) ILAY,RAD,TMAS,GS,RHO0,ISO,LAME,IPMLS
120   FORMAT(I3,F10.0,E15.7,2F10.0,3I2)
      TMASS=TMAS
C*    READ(5,130)  FMT
      READ(5,130)  FMT
C*130 FORMAT(20A4)
  130 FORMAT(A80)
      IF(ISO)  150,140,150
  140 CONTINUE
  
      DO I = 1, ILAY
      if (I.le.42) then
      READ(5,FMT)  HI(I),RHOI(I),VPI(I),VSI(I)
      else
      READ(5,111)  HI(I),RHOI(I),VPI(I),VSI(I)
      endif
      enddo
      
111   format(F6.2,F8.4, F9.4,F7.4)
      
      DO  145  I=1,ILAY
      XII(I)=0.0
      PHII(I)=0.0
      ETAI(I)=0.0
  145 CONTINUE
      GO TO  200
  150 CONTINUE
      
      READ(5,FMT)  (HI(I),RHOI(I),VPI(I),VSI(I),XII(I),PHII(I),ETAI(I),
     1              I=1,NLAY)
      
  200 CONTINUE
      READ(5,210)  NDIV,DMAX,MPUNCH                         
  210 FORMAT(I5,F10.0,I2)
C
      CALL  DIVIDE
      
      CALL  GRAVTY
      
      IF(LAME)  260,240,260
  240 CONTINUE
      
      DO  250  I=1,NLAY
      EMU(I)=RHO(I)*EMU(I)*EMU(I)
      ELAMB(I)=RHO(I)*ELAMB(I)*ELAMB(I)-2.0D+0*EMU(I)
  250 CONTINUE
C
  260 CONTINUE
      LAME1=1
      IF(JSET-7)  290,290,270
  270 CONTINUE
      REWIND  JSET
      WRITE(JSET)  NAME
      WRITE(JSET)  NLAY,RAD,TMASS,GS,RHO0,ISO,LAME1,NDIV,DMAX
      DO  280  I=1,NLAY
      WRITE(JSET)  H(I),RHO(I),ELAMB(I),EMU(I),GR(I),XI(I),PHI(I),ETA(I)
  280 CONTINUE
C
  290 CONTINUE
      
C
      IF(MPUNCH.EQ.0)  GO TO  1000
      WRITE(7,310)  NAME
  310 FORMAT(20A4)
      WRITE(7,312)  NLAY,RAD,TMAS ,GS,RHO0,ISO,LAME1
  312 FORMAT(I5,F13.5,E15.7,F10.5,F10.6,2I2)
      IF(ISO.NE.0)  GO TO  350
      WRITE(7,314)
  314 FORMAT(9H(4D20.12))
      WRITE(7,316)  (H(I),RHO(I),ELAMB(I),EMU(I),I=1,NLAY)
  316 FORMAT(4D20.12)
      GO TO  1000
C
  350 CONTINUE
      WRITE(7,354)
  354 FORMAT(17H(4D20.12/3D20.12))
      WRITE(7,356)
     1(H(I),RHO(I),ELAMB(I),EMU(I),XI(I),PHI(I),ETA(I),I=1,NLAY)
  356 FORMAT(4D20.12/3D20.12)
      GO TO  1000
C
C TAPE INPUT
C
  500 CONTINUE
      READ(ISET)  NAME
      READ(ISET)  NLAY,RAD,TMASS,GS,RHO0,ISO,LAME,NDIV,DMAX
      TMAS=TMASS
      DO  510  I=1,NLAY
      READ(ISET)  H(I),RHO(I),ELAMB(I),EMU(I),GR(I),XI(I),PHI(I),ETA(I)
  510 CONTINUE
C
      DD=H(1)
      HH=0.0
      DO  520  I=1,NLAY
      DD=DD+HH
      HH=H(I+1)
      D(I)=DD
  520 CONTINUE
      GO TO  260
C
C EXIT
C
 1000 CONTINUE
      RETURN
      END
